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COMPUTER SIMULATION OP SURFACE AND FILM PROCESSES 


In this report the results obtained in two separat 
investigations are presented: 

(i) Simulation studies based on a potential comprising two-body 
and three-body interactions. As model systems» several dif- 
ferent silica phases are taken into consideration, 

(ii) Crack propagation studies in two-dimensional triangular 
lattices , 

The progress made in these areas during this period is 
given in the following two self-contained chapters. 
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Simulation studies based on a potential with two- and 
three-body interactions. Crystalline silica phases, 


their surfaces and interfaces were modeled. 



INTRODUCTION 
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Using comptiter simulation techniques, the behavior of a model system 
(■with many degrees of freedom) can be computed to analyze various physi- 
cal phenomena. It has been shown that adequate computer methods, based 
on interactions bet'ween discrete particles, can provide information leading 
to an atomic level understanding of various processes [1-7], The success of 
these simulation methods, however, is closely related to the accuracy of the 
potential energy function representing the interactions among the particles. 
The prime objective of the present investigation is to develop a potential 
energy function for crystalline S 1 O 2 forms that can be employed in lengthy 
computer m'-deHing procedures. Such a potential energy function which 
could furnish acceptable crystalline forms of Si 02 at iBnite temper ature.s 
(say, under a Monte Carlo or molecular dynamics code) presently is un- 
available. In many of the simulation methods which deal with discrete par- 
ticles, semi-empirical two-body potentials have been employed to analyze 
energy and structure-related properties of the system [1-10]. However, it is 
now well accepted that many-body interactions are required for a proper 
representation of the total energj' for many systems. There are numerous 
reports indicating the significance of multi-body interactions in calculating 
energy and structure-related properties[ll-19j. In this study, we included 
many-body interactions in an appropriate way for simulations based on 
discrete particles. 

In general, for a system of N particles, the total potential energy may 
be expanded as [15]: 



. N N N 

' % 3 " 


( 1 ) 




where, u(ri, ry), t/(f , rjt), ..., u{ri,rj, are two, three, and n-body poten- 
tials, respectively. The position of the i’th particle is denoted by 

Clearly, the most important term in this expansion is the first term 
involving two-body interactions. Therefore, in the majority of the atomis- 
tic calculations made to date, only pairwise additive potentials have been 
used. This provides great simplification in the analytical formalism as 
well as in numerical computations. When the parameters of a two- body 
potential, however, are calculated from, say, an experimentally determined 
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bulk property, the pairwise nonadditive (Le, multi”body) contribution to 
the total energy is introduced indirectly via estimated parameters U a 
complicated way. The two-oody potential, therefore, with parameters car- 
rying the effect of the many-body portion of the total potential is called an 
"effective” pair potential [14,15,20]. In general, effective potentials proTiide 
good results when another property of the same state is calculated (i.e., 
as long as variations in the local atomic configuration remain negligible) 
[21]. However, if the propei^des of another state have to be calculated, the 
many-body part should be included in a theoretically acceptable way, if its 
contribution is considerable [20], 

In this investigation, energy and structure-related properties for sys- 
tems containing Si and O atoms were taken into consideration. Simulations 
were performed at finite temperatures using a Monte Carlo procedure. For 
total energy calculations, both two-body and three-bocy interactions were 
taken into account. Reported model calculations for euergy and structure- 
related properties of systems with Si and O atoms are rather scarce [22- , 
26]. In general, all these earlier calculations were carried out assuming 
only two- body "effective" potentials. An ionic model based on the Born- 
Meyer-Huggins potential (which is an "effective" pair potential in nature) 
was adopted for this purpose [27]. For several macroscopic properties, in 
particular for the representation of amorphous structures, this potential 
function produced very interesting results which are consistent with ex- 
periments. However, at the present time, there is no indication that this 
potential would allow various crystalline Si 02 structures to be stable [1]. 

In the present study, with a potential energy function comprising 
two and three-body interactions, we were able to obtain stable crystalline 
structures for low and high temperature forms of quartz and cristobalite 
employing a Monte. Carlo procedure. The results emphazise the importance 
of multi-body interactions in the stability of crj^stalline Si 02 phases. 

POTENTIAL ENERGY FUNCilON 

To calculate the total potential energy of a system in a given configuration, 
equation 1 was employed considering two-body and three-body interactions. 
Four and higher body terms, however, were neglected. Because this expres- 
sion (equation 1) had to be used in lengthy machine computations, ti(r,,ry) 
and.u(T,-, fj, ?*) were chosen with the simplest possible functional forms. In 
this study, therefore, the twr-body part was represented by a Mie-type 
potential which is given by: 




(m — n) 



( 2 ) 
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wlth,rij‘ = I"?! — ^j\, To denotes the equilibrium separation and e is the 
energy at ry = Te> The exponents m and n account for the repulsive and 
attractive terms, respectively. The three-body part, on the other hand, -was 
expressed as: 


“(?.■, Oif*) = Zj • c?i(r<, ry,rt) (3) 

I 

■where, the summation includes all triple multipole interactions resulting 
from the expansion of the third-order interaction energy for three atoms, 
Each term in the summation is expressed as the product of a geometrical 
factor G'(ri,fy,rfc; which depends on the relative positions of the three 
atomic nuclei and an interaction constant which depends only on the atomic 
species involved in the interaction. The functional forms of G(r», ry, r*) for 
several multiple interactions have been obtained by Bell [28]; and by Doren 
and Zucker [29]. Here, we consider only the triple-dipole interaction which 
has been shown to be the dominant contribution [30,31]. This term was 
first obtained for closed shell atoms by Axilrod and Teller [31] as: 


ry, r*) ^ Zi • ^i(ri,fy,f*) (4) 

with 




1 ZCosOi CosBj CosBi; 


( 5 ) 


where, and rijfTikiTjfi represent the angles and the sides of the 

triangle formed by the three particles i, j and k. Parameters such as e, 
m, n and Z will be referred to as ” internal” parameters indicating that they 
depend only on the atomic characteristics of the corresponding particles. 

Now, combining equations 1 through 5, and neglecting the four-body 
and higher terms, one obtains: 


$ = -r ^3 


( 6 ) 


with 
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and 
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In this general form, the total potential energy $ is a function of the atomic 
coordinates, and the internal parameters. These parameters, as mentioned 
above, are assumed to be independent of the atomic positions and the 
geometry of the system; they depend only on the atomic properties of the 
species involved. 

For a binary system of types 1 and 2 with the corresponding number 
of particles Ni and N 2 (for example, for the Si02 case), equations 7 and 
8 may be rewritten as: 


( 


1 + 3Cos$i CosOj Cos$k 
irij-rik-rjk)^ J 


( 8 ) 
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^ y 2(mafi - na^) 
with 0 = 1 , 2,, >5=1,2 and 7 = 1 , 2 ; and 
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where, 

l.’W 

R.f = ^ 

( 11 ) 
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and 





Ta fl'j = To-fc) 

J fc 

(14) 


with, Tai 7 ^ 0, r^j 0 and r^k 0; (and considering Here, 

represents' the equilibrium separation between the two species q and 


p for the two-body part ; d is a critical distance parameter which is used 
as a normalizing factor. ^ is independent of the value of d. raj = \?a — 
fjl, where, fa denotes the position of an atom of type a taken as the central 
particle for the summation. The total number of particles in the system is 
given by Ntoi == M -h N 2 ; while, rriafi and nap are the exponents for the 
repulsive and attractive parts of the two-body potential. 

The hCe potential (equation 2) which is employed to represent "bare” 
two-body interactions is the general form of the commonly used Lennard- 
Jones potential with m and n assumed to be equal to 12 and 6, respectively. 
Here, m and n determine the curvature of the function around the mini- 
mum and would provide a greater flexibility, in particular, when left as 
variables in calculating quantities expressed as derivatives of the potential 
energy. Even though the Mie potential may represent "bare" two-body 
potentials fairly well, an accurate representation of the three-body part is 
more difficult to attain. Only a limited number of reports can be found 
in the literature proposing an analytic formula for three-body interaction 
which would meet our requirements. In this first attempt, we adopted the 
Axilrod-Teller approach to analyze, at least, some of the generic features 
produced by a formal three-body function in the calculation of structural 
and energy-related properties for various Si 02 crystals. 

The functional representation of ^ as sums of two and three-body 
terms is quite different from the usual formulation of the total potential 
energy as a sum of "effective" two-body potentials alone. In the present for- 
malism u{7i,7j) is the "bare" two-body potential and the many-body effect 
is introduced via u(r,,ry, r*). In the "effective" pair potential formalism, 
however, the many-body part is incorporated with the two-body ’internal’ 
parameters, therefore, the muiti-bo/iy part was implicitly assumed to be 
pairwise additive. This assumption, in turn, makes the ’internal’ potential 
parameters very much more locally structure dependent (i.e.,for another 
geometrical arrangement of atoms, for instance, going from bulk to the 
surface, adjustments in the parameters would be necessary). At this point, 
it is anticipated that the present formalism of the total energy as sums 
of "bare" two-body plus three- body interactions would cause the internal 
parameters to become much less position and geometry dependent than in 
the case of an "effective" two-body approximation. See appendix 1. 

CALCULATION OF THE PARAMETERS 

For a binary system such as Si02, evaluation of the parameters (internal 
parameters, in particular) constituted the most cumbersome part of this in- 
vestigation. We used only experimental data to evaluate these parameters. 
The total potential energj' of a binary system as given by equations 6 and 9 



throTigh 14 contains two classes of parameters, (i) Internal parameters 
which are related implicitly to the atomic properties (such as the polarizability, 
mean excitation energy and atomic volume) of the species involved in the in- 
teraction. (ii) Geometric (external) parameters, as given by equations 12 
through, 14 are functions of atomic positions and depend solely on the 
configuration of the system. These geometric parameters are basically 
some kind of lattice sums and, for known crystalline structures, can easily 
be calculated. Table 1 tabulates these parameters for several crystalline 
S 1 O 2 forms. Internal parameters, on the other hand, are more difficult to 
evaluate. In contrast to the geometrical parameters which are dependent on 
the configuration, internal parameters are assumed to be independent of the 
habit and the geometrical setting of the system. Therefore, the iame set of 
interna! parameters which reproduces experimental quantities satisfactorly, 
say for the pure Si, can be employed for calculating the two and three- 
body interactions among Si atoms in any SiO^ form. The total potential 
energy, expressed by equation 6 is, by definition, the cohesive energy of , 
the system. This constitutes one of the main expressions in the internal 
parameters evaluation procedure. The other major relation that has been 
used is the stability relation given by: 


dV 


= 0 


(IS) 


where, V denotes the total volume of the system. Obviously, this expression 
is exact only at static equilibrium (i.e., at T =0®K). However, it has been, 
generally, assumed that the minimum of $(Y) coincides with the minimum 
of the free ener^ of the system as a function of V. 

For a binary system such as 5t02, the internal parameters can be par- 
titioned into three groups, (i) Parameters for pure silicon: Experimental 
small cluster and bidk data for Si were used to calculate pure silicon 
parameters [32,33]. Equation 6 has been employed in this evaluation as- 
suming m^si—si) = 12 and n^si—si) = 6. The two body energy well 
depth for equation 2 was calculated as e( 5 i_si) = 37740.0'’iiL based on 
the reported experimental dissociation energy of the Sj ‘2 molecule. The 
standard enthalpy of vaporization for the crystalline Si has been reported 
as 56900 i 1500®K [33]. This value was used in equation 1 as an input (a 
small term of the order of PV has been neglected, here, and cohesive energy 
has been assumed). Then, the solution of equation 1 based on 
the stability criterion (equation 15) produced: To{Si—si) = 2.2134 A and 

Z{si—si—si) = 5.2934 X 10’^(®iiCA ). This set of parameters reproduces 
the correct energy and structure for the crj^stalline silicon, as well as the 
reported m'icrocluster data (i.e., energy and the structure) for 5*2 and Si^. 
Furthermore, the diamond structure for the bulk Si satisfies the stability 



criterion and was found to have a lower energy than either the hep, fee or 
bee structures. These parameters, however, were not tested for reproduc- 
ing quantities expressed as higher deri-^utives of (ii)Parameters for pure 
oxygen: Considering the dissociation energy of the O2 molecule, the two- 
body potential well depth for equation 2 was taken as f(o— o) = 59943.0 
(^K) [34]. The exponents in the two-body function were assumed to be 
m(o— o) = 12 and n^o—o) == 6- Bas^ml on spectroscopic data the equi- 
librium distance between two O atoms was taken as the average bond 
len^h for the O 2 molecule, rr(o—o) == 1.208 A [34]. In this case no bulk 
data (for example, based on solid oxygen) was employed to calculate the 
three-body intensity parameter. Instead, we used reported ozone (O3) data 

[35] and obtained 2’(o— o— o) = 2.622 x 10® ("iCA*). These parameters 
for pure oxygen are, now, exact for the minimum point of the two-body 
potential curve and provide good dissociation energy for O3. However, 
minimization of the total energy yields a linear O3 molecule. At this point 
no further attempt has been made to correct this shortcoming by refining 
the parameters, (iii) Cross parameters; in the calculation of the cross 
parameters (i.e., for interaction between Si and O atoms) experimental 
data from bulk o-cristobalite and a-quartz were employed. In this case, 
again by assuming m{si—o) = 12 and n^si^-^o) == 6, we are left with four 
cross parameters (€(^^si^o),ro{si^o), ^(si-si^o) and Z(si-o-o) ) to be 
determined. The cohesive energies for a-cristobalite and a-quartz were ob- 
tained from standard enthalpies of formation which are tabulated in Table 
2. These experimental cohesive ener,gies for a-cristobalite and a —quartz 
were employed in equation 6 in determimng the parameters. Now, equa- 
tion 6 and 15 for a-cristobalite and a-quartz constitute four "nonlinear’’ 
simultaneous equations. However, in these simultaneous equations, the 
coefficients for the corresponding terms are numerically very similar, and 
this makes acceptable roots difficult to obtain. To overceme this undesired 
"ill condition" situation and, at the same time, to investigate the effect 
of the input values of cohesive energies on the parameters, the 'V'alues of 
4^’s for both a-cristobalite and a-quartz were varied, within limits of the 
experimental error margin. It was found that, in general, small variations 
in ^ produces acceptable roots. So far, we have considered only positive 
roots. Table 3 tabulates calculated values for the cross parameters using 
two different sets of cohesive energies (varied within experimental limits), 
along with the parameters for pure s3^stems. 



COMPUTATIONAL PROCEDURE AND RESULTS 

Emplojing parameters obtained in the preceding section, simulation 
calculations were carried out at finite temperatures for cv^crlstobalite, o- 
quartz, ^-cristobalite and ^-quartz. (These four allotropic crystalline forms 
of Si 02 will be referred c and P — q, respectively). 

In this investigation we are mostly concerned about the structural stability 
of these four Si02 forms under a Monte Carlo code based on a potential 
function given by equation 6. The stability criterion imposed by equation 
15 in calculating internal parameters may only be valid at the static limit 
(i.e., at T= O^K), This condition, therefore, can not necessarly guaranty 
the preservation of the proper symmetry and configuration of the system 
at a finite temperature. 

Model calculations were performed using a Monte Carlo technique 
based on the Metropolis approximation [36-38]. This is a stochastic method 
weighted according to the Boltzmann factor. In this procedure, every 
particle is treated discretely at finite temperatures. In general, the particles 
were positioned initially according to the crystalline structure of the system 
under consideration. Then, every atom was displaced by a self-adjusting 
step size. The acceptance of the new position is decided by the Metropolis 
procedure which is employed for every particle in the system based on the 
temperature and the particle’s total energy. This technique generates a 
series of snapshots representative of the energetically most probable region 
of the phase space and, in principle, it can provide ensemble averages of 
any position-dependent quantity. 

Throughout this investigation, calculations were performed consider- 
ing the same internal paramet-ers for pure silicon and pure oxygen. In the 
case of the cross parameters we used set #1, as given in 'Table 2; occa- 
sionally, however, set #2 was employed for parallel runs to investigate the 
effect caused by the cross parameters on the final result. This analysis, 
namely, the investigation of the effect of the parameters on the calculated 
macroscopic characteristics of the system, at the moment, is incomplete. 
However, the present study indicates that they are extremely important in 
understanding the intricate relationship between atomic properties and the 
structural characteristics of the system. 

First, the bulk structures of the four different Si02 forms (a— c, cr — g, 
^ — c and P—q) were simulated. The temperature was taken as 300 ’K and 
the initial configuration for each system was generated based on its reported 
crystallographic data [39-41]. For calculating bulk properties, a three- 
dimensional periodic boundary condition was imposed on the systems to 
eliminate the surface effect [42]. The volume of the computational box was 

• 3 

taken somewhait larger than lODOA depending on the crystalline structure 



of the system. The potential energy cut-off radius, Rcuti 'ffas taken typically 
as 5.0^. In order to analyze the effect of the on the final result, the 
calculation for o — c i>ras r^ipeated with Rcut — 7.0A, but no appreciable 
difference in either the total energj' or in the final atomic configuration 
was detected. All iterations were carried out until complete equilibration, 
which was monitored by the variation in the total potential energy. In 
general, equilibrations were attained before 5000 iteration steps. For all 
the four SiO^ forms studied, it has been found that equilibrated structures 
preserved their original symmetry very well. The orientation of atoms 
in these equilibrated structures was analyzed by the radial distribution 
function (RDF), and by the angular distribution function (ADF). 

The RDF provides mean peak positions for d(^si—o), rf(o— o) and 
d(5i— S i) corresponding to 'Irst neighbor distances between Si— 0, 0 — 0 
and Si — Si atoms. Ai an example, figure 1 demonstrates the calculated 
total RDF for the a — • c case. At 300® A, the mean peak positions for all 
four cases remained virtually unaffected (with respect to the initial first 
neighbor distances' based on the crystallograpliic data). Due to thermal 
fluctuations, however, the peaks became more diffuse, as anticipated. The 
most pronounced relaxation took place for the Si — 0 peak of ^ — c for 
which we considered the ideal structure as the initial configuration. The 
equilibration, in this case, caused a shift for d(c,-._o) fr^ni 1,71 to 1.57 (to 
1.58 for the cross parameter set #2 ), which is in, the right direction toward 
the experimental value of 1.612 for the real — c [39]. Furthermore, the 
RDF mean peak positions were found to be quite insensitive to the type of 
cross parameters employed in the relaxation procedure (see Table 4). 

The ADF was used to analyze average Si — 0 — Si and O — Si — O 
angles formed by the nearest neighbors of the corresponding atoms. These 
angles are found to be very useful in the characterization of various 5:02 
forms [25,39]. In all four cases investigated here, averaged values of the 
O— Si— 0 angle, which is the internal angle of the Si 04 tef rahedron, was 
found to be around 109 degrees. This is consistent with many experimental 
findings. The intra- tetrahedral angle, Si — 0— Si, however exhibits some 
’V’ariations among the allotropic forms of Si02 considered. As an example. 
Figure 2 demonstrates the ADF for Si — O—Si and O •— 5i — O of a — c. 
The effect of the type of cross parameters was found to be somewhat more 
visible on the average value of the Si—O — Si angle when compared with 
the effect on the averaged 0—Si— 0 angle. In general, cross parameter set 
#1 provides slightly larger Si—O—Si angles (see Table 4 for comparison). 

Calculated cohesive energies for the equilibrated 5:0*2 forms are given 
in Table 2, '^long with the experimental values. All the calculated cohesive 
energies were found to be within acceptable limits. The effect of different 


cross parameters on the calculated total energies is noticeable but quite 
small and the calculated energies remained vithin esperimental error mar- 
gins. 

Next, we performed simulation calculations for (001) surfaces of the 
a — c. Again, the Monte Carlo method based on equation 6 was employed 
considering the same set of parameters, as above. Initially, exposed surfaces 
were generated by cutting an, ideal a— c crystal in the (001) direction. This 
was accomplished by imposing periodic boundary conditions only in the x 
and y-directions for a properly positioned crystaJ that produced "infinite^ 
exposed surface planes in both and — z directions. One obtains a 
"regular” surface structure when the cutting plane is located between two 
neighboring 0-planes in the (001) direction. Otherwise, depending on the 
location of the cutting plane, the surface created would be O-rich in one 
side and 5i-rich in the other. Equilibration runs were carried out (i)for a 
regular surface and (ii) for an O-rich surface, considering systems with 96 
moving atoms at T= 300®iC. 

The reconstruction for the "regular" surface was found to be minimal. 
The largest displacement was^exhibited by the top O atoms which expanded 
outwards approximately 0.7A At the same time these top layer 0 atoms 
moved slightly toward the nearest Si atoms. Top and side views for 
the regular' ?>mface are shown in Figure 3 before and after equilibration. 
Equiilbrs^ted surface configurations obtained using cross set #l and set 
#2 were virtually identical. However, the effect of the cross parameters was 
very pronounced in the surface energy calculations. Cross parameter set 
#1 produced a surface energy of 841 .6 mJlm^ which is in a good agreement 
with experimental findings [43], Set #2, on the other hand, yielded 
a negative surface energy value (—310.0 mJJm^) which is in obvious 
disagreement. 

’ * The O-rich surface for the a— c exhibited considerable surf ace reconstruc- 
tion. O atoms located in the top layer (i.e,, excess O atoms) moved closer to 
0 atoms of the neighboring tetrahedron (located slightly lower than the ex- 
cess 0 atoms) forming so called peroxide bridges. Figure 4 shows the 
schematic representation of the reconstructed O-rich (001) surface of a — c. 
Formation of these peroxide bridges has already been determined by several 
investigators [44,45]. 

DISCUSSION AND CONCLUSIONS 

The effect introduced by three-body interactions is quite pronounced 
and complex. In contrast to two-body interactions, three-body forces en- 
courage open structures. Therefore, in the present formalism the familiar 
picture of pairwise additive interacting particles, which favors compact 


structures; is drastically altered. For example, the behatlor of the 0 atoms 
in the Si02 crystals represent this three-body eflfect quite clearly. The in- 
ternal parameters used for the pure oxygen two-body potential were chosen 
to reproduce the equilibrium O 2 data correctly. Accordingly, in the two- 
body potential part for pure 0 we have ro(o~o) = 1*208 A (see above). In 
the simulation runs for crystalline Si 02 forms, the average nearest neigh- 
bor O -- O distance was found to be around 2.59 A . Strong three-body 
interactions, due to surrounding Si or other O atoms in the bulk prevented 
the O atoms from further approaching each other (despite thermal motions 
introduced by the Monte Carlo code employed). At the 0-rich surface, on 
the other hand, due to a "reduced" three- body effect, the top 0 atoms can 
come closer to form peroxide bridges. 

The present interaction potential function does not utilize any covalency 
or ionicity concepts, nor does it catagorize the interactions between particles 
into short or long range. This concept, although it may be quite useful in 
many other modelling procedures, would be rather difficult to assess, in par- 
ticular, at the surface region oi in the vicinity of a defect. 

One of the most important outcomes of the present investigation is the 
fact that three-body potentials, even in a simple analytic form (but, ade- 
quately represented by triplet summation), can provide stable Si 02 crystal- 
line forms involving thermal fluctuations. This work, in this respect, may bn 
regarded as a first attempt at investigating the multi-body effect in the ener- 
getics and structure-related properties of St02 phases. Calculated potential 
energy parameters reported in this paper represent only a sample case for 
comparison. The test for the ability of these parameters to reproduce other 
macroscopic quantities of the corresponding compounds is an enormous job 
that was only performed partly in this investigation. Therefore, reported 
numerical results may be only semi-qualitative in nature, and should be 
treated accordingly. 

The computational time required for these simulation runs were con- 
siderably higher than the usual calculations based on two-body interac- 
tions alone. For example, an equilibration run for a system containing 
approximately 200 relaxing particles may take up to an hour of CPU time 
in a CDC 7600 machine. The largest part of this computation time is 
spent to calculate the three-body interactions. Despite this relatively long 
computer time requirement, simulation runs for Si 02 systems based on the 
present potential energy function are well within reach. 



APPENDIX 1. 


In tho present model the ’bare’ tTVo-body potential, for a par- 
ticular pair is always represented with the same functional form associated 
with the same parameters irrespectiTe of the immediate surrounding. The 
effect of the neighboring particles (which is the many- body effect) is ac- 
counted for via the three-body term. Therefore, the evaluation of the 
two-body part can be accomplished using simply diatomic data. However, 
extreme care must be exercised in extracting such information from ap- 
proximate ab initio calculation results. 

For any given two atoms the "bare” two-body potential is a function 
of the internuclear separation, nj, and formally it may be defined as: 

~ ^iUjoo) Al 

where, E(r,-y) and E[rijoo) denote the total ground state energies for the 
two particles at a separation Uj and at an infinitely large internuclear 
separation, respectively. Based on first principle quantum mechanical con- ‘ 
siderations, E(r,-y) and E{rijoo) may be exactly calculated by: 

E{Ui) = Al 


B(r,;oc) = AS 

■where, and are the total Hartree-Fock and electronic 

correlation energies for the two particles at a separation r,-y, and 
E^corr)^^^'^ and E^^^\2), denote corresponding energies for 

the isolat-ed particles 1 and 2, respectively. In the absence of any ex- 
ternal field, both E{rij) and E{rijoo) represent the lowest energy level 
for the corresponding systems. Accordingly, the function i/(rty) can be 
derived continuously from the long range dispersion limit to the short 
range steep repulsion region which also covers the region around the min- 
imum. Obviously, the quantum theory for such a curve does not re- 
quire separate approaches for the long and short range regions. 

In particular, for many electron systems, an exact evaluation of w(r,j) 
based on first principles, at the moment, is not possible due to various 
computational difficulties. To overcome some of these inherent difficulties, 
various approximate methods have been proposed. In the majority of 
these methods for calculating u(rij) involving many electron systems, one 
generally assumes a particular hybrid type for the bonding between the 
atoms. Accordingly, calculations produce a u(r,y, /?) curve which is a two- 
body potential curve for the particular hybrid type imposed. For different 


' hybrid types one obtains a series of u(r,y, h) curves some of vrhich may 
coincide "with w(r,'y) for a certain range of r,y. In principle, the low 
envelope the these u(r,y, h) curves must coincide vith u{t{j) of equation 

Al. 
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FIGURE CAPTIONS 

FIGURE 1. Radial Distribution Function for o-cristobaiite at 300®iC. 
The curve is obtained from the positions of Si and 0 atoms averaged over 
500 steps after full equilibration. 

FIGURE 2. Distribution Function of 0 -- 5t — O and — O — Si 
angles for a-cristobalite at 300® iiC. The results represent averaged values 
over 500 steps after full equilibration. 

FIGURE 3. Schematic representation of the regular (100) surface of 
a-cristobalite; (a) before and (b) after relaxation. The top view shows 
one of the periodic cells from the z-direction (perpendicxilar to the exposed 
surface). The side view shows the same cell viewed from the x-direction. 
During the relaxation process, only the top 0 atom was displaced slightly, 
while the rest of the system remained virtually unchanged. (Large open 
circles represent Si atoms, and small solid circles are O atoms,) 

riGURE 4. Schematic representation of the reconstructed 0-rich (001) 
surface of a-cristobalite. Peroxide bridges formed by excess 0 atoms are 
responsible for the surface reconstruction. {Si and 0 atoms are shown in 
large open and small solid circles, respectively.) 
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CHAPTER II. 


Crack propagation studies in two-dimensional 
triangular lattices. 


PROPAGATION OF CRACKS IN 2D LATTICES 


This investigation is an extension of the earlier study which 
had been carried out last year considering a static approach. 
General characteristics of the model considered in this new study 
are basically similar to the previous one. Again a 2D triangular 
lattice has been used with particles interacting via Lennard-Jones 
type potentials. However, in this case, the effect of the temper- 
ature was taken into account by using a molecular dynamics simu- 
lation technique. Therefore, the results are expected to be more 
realistic. Simulations for two-dimensional systems are rela- 
tively easier to analyze than results for three-dimensional 
systems. First of all, 2D systems contain a smaller number of 
particles and, therefore, require less computer time. Results can 
be represented by simple 2D plots and problems arising due to the 
multi-particle character of the system are easily identifiable. 
Thus, the analysis of the 2D system provides considerable infor- 
mation not only about the microscopic nature of the crack growth 
phenomenon, but also provides some knowledge about "how to 
interpret the simulation results." The question of the 
credibility of these 2D results, of course, remains unanswered. 

At this stage, it is not known how to extrapolate results obtained 
from a 2D simulation to a 3D domain. However, the results obtained 
in this study together with several other reports [1,2] in the 

i 

literature indicate that 2D systems, in most cases, do exhibit 
characteristics similar to 3D systems. 

Atomistic level ana:,|’sis of the crack propagation process 


using conputer simulation techniques has been the subject of 
several earlier investigations. In the literature we could find 
only a few reports relevant to the study carried out in this 
investigation. In the report by Ashurst and Hoover [1], the 
fracture phenomenon was investigated based on a truncated Hook's 
law force. They have found that, even with this very simplistic 
force law» their static simulation furnished results for energy, 
entropy, stress concentration and crack structure all to be con- 
sistent with expectations from macroscopic elastic theory. 

The other relevant and more recent study was reported by 
Dienes and Paskin [2}. In this modeling study they also consi- 
dered a 2D triangular lattice with particles interacting via the 
Lennard-Jones function. A crack has been introduced in the 
interior of a pre-stressed sample. The crack was initiated by 
"cutting" the bonds between a given number of atoms at the central 
portion of the sample. The interatomic potential was artificially 
set to zero between these atoms. According to their report, the 
condition would correspond to the insertion of a very thin knife 
to create the crack. Furthermore, in the energy and force calcu- 
lations, they only considered nearest neighbor interactions (by 
taking Rcut = 1.6 ro). In their model, the crack was aligned 
parallel to close-packed rows and displayed a linear path in its 
propagation. Finally, they found that their results are quantita- 
tively good at the early stages of the propagation process. 

The main objective of the present study is to investigate 
the crack propagation phenomenon at a microscopic level. For this 
purpose, a 2D triangular lattice was taken into consideration and 
the effect of a tensile load imposed on the system was analyzed. 



In the following sections, first a perfect lattice, then a lattice 
with a surface crack under various conditions was investigated. 


PERFECT LATTICE 

As a perfect 2D lattice, the basal plane of an hep lattice 
was taken into consideration. A system of 2J^00 particles in a 
rectangular shape (80 by 30) was first generated in static 
equilibrium. A tensile load was applied in the [112] direction, 
which is the close-packed direction. This direction was also 
chosen as the x-direction in our cartesian coordinate system. 

The load was imposed in small Incremental strains (in this case 
elongations) of 0.01. This was performed uniformly throughout the 
system by factorizing all the x-compononts of the position vectors 
describing the system. In the x-direction, periodic boundary 
conditions (PBC) were applied to provide continuity for the system 
(in the tensile direction), and also to furnish two free surfaces 
in the y-directior,.. In a general sense, the imposed PBC provides 
the desired tensile strain on the system. 

The system was relaxed after every incremental strain by a 
molecular dynamics code. A cut-off radius. Rcut, of 2.86 ro was 
considered for the energy and force calculations. This Rcut is 
between the fourth and fifth neighboring shells surrounding the 
central atom and provides 30 neighbors. The reduced time step was 
taken as 0.01 and the reduced temperature was T« = 0.02 (to com- 
pare with real systems; e.g., for copper these represent 2..5E-15 
seconds and 100 K, respectively. Additional information about the 
cut-off radius and the molecular dynamics program are included in 
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Appendix I and II. The stress-strain curve calculated as a result 
of this elongation-relaxation process is shown in Fig. 1 up to e = 
0 . 88 . 

2D LATTICE WITH AN EXISTING CRACK 

A lattice with an initial surface crack was generated by 
removing 9 particles from the lower surface region of a perfect 
lattice (see Fig. 2). This system^ now with an existing surface 
crack, was elongated and relaxed by the molecular dynamics code in 
a similar way explained above for the perfect case. First, the 
effect of the temperature on the stress- strain curve was analyzed 
up to e = 0.03. Figure 3 shows two curves, dotted and solid, 
representing the stress-strain curves for T* = 0.1 and T* = 0.02, 
respectively. The shift in the dotted curve (high temperature 
curve) is mainly due to the thermal expansion. For lower strain 
values, these curves represent fully equilibrated systems. 

However, for strains higher than 0.02, systems may require 
additional relaxation times to equilibrate completely. The 
difficulty involved in attaining the equilibrium is mainly due to 
large fluctuations displayed by the stress values calculated as 
derivatives of the energy. At this stage, we believe that the 
general trend exhibited by these two curves is sufficiently 
accurate for the present investigation. Any further incremental 
elongations (in addition to e s 0.03) cause the crack to 
propagate. Determination of the critical strain; i.e., the strain 
at which the crack first starts to propagate, is difficult to 
assess. For this purpose we performed three separate runs with 
three different pre-stralned systems, namely with e = 0.03, 


e = 0.035 and e = 0,Qi\, all at T* = 0.02. The 2D lattice with the 
surface crack was strained in one single step from its original 
length up to 3.0» 3.5 and 4.0^ elongations. In the case of e s 
0.03» the crack did not exhibit any growth and the overall 
configuration of the systems remained unchanged up to 3500 time 
steps. However» for both e = 0.035 and 0.04 cases, the 
propagation of the crack took place. In these prestrained cases, 
we simulated the system under nonequilibrium isothermal 
conditions. For the e = 0.035 casi, the crack growth first 
initiated after 1000 iteration steps. Figure 4 displays the 
stages of this relaxation process up to 3200 iterations, at which 
the system reached almost to an equilibrium state. The darker 
circles in the figures represent particles with higher stresses. 
For the e= 0.04 case, on the other hand, the crack propagated much 
earlier (obviously because of the high strain imposed initially). 
The crack started growing first at the 500'th iteration step and 
the system attained an equilibrium state at approximately 2400 
iteration. The stages of this propagation process are shown in 
Fig. 5. AgUin, the darker circles display particles with higher 
stresses. In both cases, the particles at the crack tip exhibited 
high stresses consistently. Furthermore, the crack propagated 
along the close-packed rows of the lattice and, at the same time, 
tried to remain perpendicular to the applied load by choosing a 
zig Ze\g path. These behaviors are very much consistent with 
experiments and theories based on macroscopic considerations and, 
therefore, indicate the adequacy of the present atomistic level 
simulation procedure. The relaxation of the system can be 



followed in Fig. o where the average stress is plotted versus the 

V 

iteration steps. The oscillatory behavior of this curve is a 
temperature effect mainly due to vibrational motions displayed 
by Individual particles in the system. From the snapshots shown 
in Fig. 5» we also calculated the velocity of the crack propaga- 
tion. The curve in Fig. 7 represents the variation in the crack 
propagation velocity as a function of the calculated average 
stress. The upper range of this curve is near the velocity of 
sound propagation. This is expected according to a report by 
Ashurst and Hoover [1]. 

REFERENCES 

1. W. T. Ashurst and W. G. Hoover» Phys. Rev. B JJL, ( 1976 ) 

2. G. J. Dienes and A. Paskini in ”Atomlstics of Fracture," ed; 
R. M. Latanision and J. R. Pickens (Plenum Press, New York, 
1983 ) p 671 . 

3. P. 0. Esbjorn and E. J,i Jensen, J. Phys. Chem. Solids 3 7 . 
1081 (1976). 





oRieir^AL p^m ^3 

OF POOR QUALITY 

APPENDIX I 

Many of the characteristics of a 2D Lennard-Jones lattice are known 
[3], Here, the energy and its derivatives of a 2D triangular lattice will be 
formulated. For a perfect lattice, with particles interacting via two-body 
central potentials, the total interaction energy of the bulk system can be 
expressed as; 


♦ = yE“M W 

where, N is the total number of particles in the system, and u(r j) represents 
the pair potential function between two particles in the system separated 
by a distance n> For the Lennard-Jones case we have; 

«(r<) = e((^P-(^)*) (2) 

with € and to denoting the energy and the internuclear distance at the 
equilibrium, respectively. The stability criterion for a 2D system, in the 
absence of an external force, is given by; 


dA 


= 0 


where, A is the total area occupied by the system, 
particles we have: 


(3) 


For a system of N 


A = iV’a. (4) 

with Op d enoting the area per particle. For a triangular lattice is equal 
to y/o.75d^ where d denotes the nearest neighbor distance. Considering 
lattice sum notations, from equations 1 and 2 we obtain: 




(5) 

where, 



• • 

bi2=i:(V 
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(6o) 

and 

B. - 
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(64) 
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[3]. Here, the energy and its derivatives of a 2D trianguiar iattice wiii be 
formulated. For a perfect lattice, with particles interacting via two-body 
central potentials, the total Interaction energy of the bulk system can be 
expressed as; 
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where, N is the total number of particles in the system, and ti(fi) represents 
the pair potential function between t^c particles in the system separated 
by a distance r«. For the Lennard-Jones case we have; 

= ( 2 ) 

with € and r« denoting the energy and the internuciear distance at the 
equilibrium, respectively. The stability criterion for a 2D system, in the 
absence of an external force, is given by; 


dd 


= 0 


where, A is the total area occupied by the system, 
particles we have: 


( 3 ) 


For a system of N 


A = Na. (4) 

with a« denoting the area per particle. For a triangular lattice is equal 
to ^0.75d^ where d denotes the nearest neighbor distance. Considering 
iattice sum notations, from equations 1 and 2 we obtain: 

where. 


and 

^ A 

= Dr)‘ 



( 66 ) 
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are the 2D lattice sums. Now, the condition for the minimum energy 
expressed by equation 3 gives: 

7“ I'l 

The definition of the equilibrium condition based on this equation is very 
important for calculating the nearest neihgbor distance d with respect to 
the equilibrium distance rp of equation 2. The value of fo/d depends on 
the total number of particles included in che summations of equation 6. A 
consideration of a finite cut-off radius, Rcut, in the lattice sum calculations, 
therefore, would affect the values of rp/d as well as $ Accordingly, Rcut 
becomes a part of the potential function, in particular, for smaller cut-off 
radii. 

In the computer simulation studies the stress and the elastic constants 
are calculated at the atomistic level for every individual particle. The 
Lagrange strain parameters, Cap, for a 2D lattice may be expressed as; 


ff = 2? -f 2c,xJT? -f 2cyyy? -f 2cxyXiyi -f 2cy*j/,-2< (8) 

The stress, for example, in the x-direction is given by: 

_ 1 

“ ak: 

or, 




=«ir 

ft 




Op "I n 

Similarly, for the y-direction we have: 
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( 9 ) 


( 10 ) 


The elastic constants, of a 2D triangular lattice have already been 

calculated by Esbjorn and Jensen [3], and following relations have been 
obtained: 


and 


Oatxx — ^yyyy — 


xxyy 


(lla) 


Oxxxy — ^yyyx — 0 . 


where, 


(116) 



As it was stated above for r«/(/ and ^ cases the values of Cotfi and Cop^e 
are also dependent on the cutroif radius considered in the summations. Of 
course, for a system at the equilibrium (i.e., in the absence of external 
forces and if the condition described by equation 3 is satisfied) the stress 
components tap should be equal to zero. In investigating the effect of Bcuh 
it is very helpful to analyze the system in terms of neighboring shells. Table 
1 gives these neighboring shell distances and the occupation numbers. Table 
2, on the other hand, tabulates calculated values for the lattice sums, r«/d, 
^ and Cssxx upto 20’th shell. 



APPENDIX IL 


COMPUTER SOFTWARE (TENST) 

This package has been prepared mainly for investigating behaviors of 
tYio or three-dimensional monatomic systems under tensile loads. Simulations 
are based on a molecular dynamics technique -which uses a Nordsieck- 
Gear algorithm. The program -was prepared in a vectorized form for the 
Cray usage. It can handle sytems containing up to iOOOO particles within 
reasonable computation time. The result of the calculations for positions 
and individual stresses can be stored for any desired intervals for further 
analysis. This program has also a graphics option for trajectory plots. If 
desired trajectories of every moving particle can be plotted for a visual in- 
spection. 

INPUT DATA AND THE DEFINITION OF PARAMETERS 

(i) I/O data: 

1) TITLE (Title card up to 80 characters) 

2) XTIME,NRUN,NCO.NT,NEW,rPSTEP,XMlN,XMAX 

3) IFREQ,IMEAN,ISTRES 

4) NSTEP,NMEAN,IPDATE,NSBSV;NSBTS,NFINTS 

5) XELO,NELO,JDEVl 

6) TEMP,DT,RCUT 

7) EP,RO,AMASS, MOVE, FAC, LAYER 


X repeat LAYER times 

8) W(J,K),NOO(J,K),ISQUAR 


9) PP(1), PP(2), PP(3), FACPBC 
(ii) Definition of parameters: 

TriLE= Title card up to 80 characters. ' 

XTIME= Time limit for internal check. It should be less than the 
time allocated for the program run for normal termination. 

NRUN= Maximum number of steps for the run. 

NCONT= 1 For initial runs. 

NCONT= 2 For continuation runs. 

NEW=.0 

IPSTEP= Trajectory plotting routine will be called in every IPSTEP 
times. If IPSTEP=0 no plot calls will be made. 



XNIIN; XMAX iMinimum and maximum points on the plot (even if 
PSTEP^O XMIN and XMAX must be defined). 

IFREQ, INIEAN : In every PREQ steps, averages of the individual 
stresses will be calculated for the last MEAX steps. Always IMBAN.LE.PREQ 
should be satisfied. 

ISTRES= 1 Averages for stresses will be taken. 

ISTRES= 0 No stress calculations will be performed. 

NSTEP = Step number at the start (usually NSTEP=0). It has no 
effect if NCONT==2. 

NMEAN= Number of steps averages to be taken. This is for the 
periodic printnig only. 

IPpATE= Number of steps between neighbor list update. 

NSBSV= Number of steps between savings of variables (coordinates • 
and averaged values). 

NSBTS= Number of steps between temperature scalings. 

NFINTS= Number of step at which temperature scaling will be turned 
off. 

Elongation factor to be multiplied by the x component of 
the positions for incremental elongation. (If XELO=l no elongation will 
be imposed). 

NELO= Number of steps between elongations 

JDM= Dimensionality of the system. (For a two-dimensional system 
JDIM==2). 

TEkip= Temperature (in deg. Kelvin). 

DT= Time steps in seconds. 

RCUT= Cut-off radius (in A) 

EP= Energy parameter, epsilon (in deg.) 

. RO= Equilibrium distance for two-body interaction (in A). 

AMASS= Atomic mass (in atom gram). 

MOVE= Number of mobile particles. (In general MOVE should be 
equal to the total number of particles). 

FAC= Factor used in generating lattice points. 

LAYER= Number of input files required for lattice generation. 

W(J,K)= First three describe the coordinates of the original atom. 
Following three are the unit cell dimensions. 

NOO(J,K)— Number of particles to be generated in x,y and z direc- 
tions. 

ISQUAR= 2 Takes square root of the W(J,K) values before the lattice 
is generated. If ISQUAR.NE.2 no action will be taken. • 

PP(1)= Length of the periodic box in the x-direction. This is for the 
PBC option and must be consistent with the genereted lattice coordinates. 

PP(2)= 0. 



PP(3)= 0. 

FACPBC= Factor for PP(1) to provide the correct length for PBC. 


TABLE CAPTIONS 


S 


1. Neighboring shells# number of particles in every shell# total 
number of particlesi the distance squared and distances of 
shells from a central particl:-; are given for a two-dimensional 
triangular lattice, 

2. For every shell the values of lattice sums (Bg and B^g» 

total potential energy $ and the elastic constant are 

tabulated for a two-dimensional lattice (see Appendix I). 
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FIGURE CAPTIONS 


1. Stress-strain curve for a perfect 2D triangular lattice. 
Stress values indicate average stresses calculated for the 
system (in reduced units). 

2. A system with a preexisting surface crack at equilibrium. 

3. Temperature effect on the stress-strain curve for a system 
containing a surface crack. Solid and dotted lines represent 
T* = 0.02 and T» = 0.1 cases, respectively. 

il. Stages of the crack propagation for a system with a surface 
crack pre-strained for e = 0.035. Darker circles represent 
particles with higher stresses. Numbers at the right side of 
snapshots are the iteration steps. 

5. Stages of the crack propagation for e = O.O^I. See caption 
for Fig. 

6. Variation in the average stress as a function of the iter- 
ation steps. For the e = 0,04 case. 

7. Crack propagation velocity versus the average stress. (Units 
for the velocity are calculated for copper.) 
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